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ABSTRACT 

If binary intermediate-mass black holes (IMBHs; with masses between 100 and IO'^Mq) form in dense stel- 
lar clusters, their inspiral will be detectable with the planned Laser Interferometer Space Antenna (LISA) out 
to several Gpc. Here we present a study of the dynamical evolution of such binaries using a combination of 
direct A^-body techniques (when the binaries are well separated) and three-body relativistic scattering experi- 
ments (when the binaries are tight enough that interactions with stars occur one at a time). We find that for 
reasonable IMBH masses there is only a mild effect on the structure of the surrounding cluster even though the 
binary binding energy can exceed the binding energy of the cluster. We demonstrate that, contrary to standard 
assumptions, the eccentricity in the LISA band can be in some cases as large as 0.2-0.3 and that it induces 
a measurable phase difference from circular binaries in the last year before merger. We also show that, even 
though energy input from the binary decreases the density of the core and slows down interactions, the total 
time to coalescence is short enough (typically less than a hundred million years) that such mergers will be 
unique snapshots of clustered star formation. 

Subject headings: Black hole physics, gravitational waves, stellar dynamics, methods; N-body simulations 



1. INTRODUCTION 

The existence of intermediate-mass black holes (IMBHs; 
masses M ^ 10^"^ Mq) is not as certain as that of stellar-mass 
or supermassive black holes because there is as yet no con- 
clusively established dynamical mass for any candidate, al- 
though there is strong circ umstantial evidence for this mass 
range in several cases (see lMiller & Colberj|2004l and refer- 
ences therein for a review). Mergers of IMBHs would, how- 
ever, be strong sources of gravitational waves. 

The best studied scenario is the runaway growth of a star in 
a young cluster via physical collisions among the most mas- 
sive stars in the center, which have sunk through mass seg- 
regation jPortegies Zwart & McMillanI 120001: iGiirkan et al.l 
I2004t iPortegies Zwart et al."2 004t iFreitag et al.ll2006h . Re- 
cently, Giirkan et al. (2006) addressed the same configuration 
but added a fraction of primordial binaries to the stellar sys- 
tem. Using a Monte-Carlo stellar-dynamics code, they found 
that not one but two very massive stars grow in rich clusters 
in which 10% or more of stars are in primordial hard binaries, 
suggesting the formation of two IMBHs. However, this re- 
sult has not been con firmed yet using more accura te direct 
A^-body simulations. IPortegies Zwart et al.l ( |2004() have a 
simulation with primordial binaries but they do not see this 
formation, though it is also currently unclear how different 
core concentrations will affect binary IMBH formation with 
a certain fraction of primordial binaries. It is also possible 
that wind losses may drive away mass more rapidly than it 
accretes thr ough further collisions (see Belkus. van Bever, & 
Vanbeveren 120071) . although this relies on uncertain extrapo- 
lations from the ~ 120 Mq that is the top of their range (see 
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their Table 2) to the 
simulations. 
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Fregeau et al. ( |2006|) considered for the first time the pos- 
sibility that such a binary could be observed thanks to the 
emission of gravitational waves in the coalescence phase and 
estimated that one can expect the Laser Interferometer Space 
Antenna (LISA) to detect tens of them depending on the distri- 
bution of cluster mas ses and densities. Amaro-Seoane & Fre- 
itag (12006) addressed the evolution of a binary of two IMBHs 
formed as the result of the collision of two independent stel- 
lar clusters and followed the parameters of the binary orbit 
down to the region in which it will emit gravitational waves 
in the ~ 10"'^- 10"' Hz LISA domain. To do this, they com- 
bined direct-summation simulations with an analytical model 
to evolve the binary from a point in which it was hard. 

Here we assume that an IMBH binary has been produced in 
a single dense stellar cluster, and study the subsequent sinking 
of the IMBHs and the evolution and properties of the binary 
when it forms. In § 2 we discuss our numerical method, which 
combines direct A^-body studies with three-body scattering in- 
tegrations. In § 3 we discuss the astrophysical implications of 
our results. 

2. EVOLUTION OF THE IMBH PAIR: NUMERICAL METHOD 

2.1. Direct N—body simulations 

Direct A^-body codes integrate all gravitational accelera- 
tions in a stellar system without supposing any special sym- 
metries. They are thus the most gener al and robust tools fo r 
numerical analysis of stellar clusters (lAarsethlll999L l2003l) . 
The code we use, NBODY4, includes a variety of sophisti- 
cated approaches that improve speed and accuracy, including 
KS regularization ( Kustaanheimo & Stiefel 1965), as well as 
triple (3-body subsys tems), quad (4-body subsystems), and 
chain regularization (lAarsethlll999L l2003l) . It also does not 
make use of any softening, which would lead to unrealistic 
evolution of the orbital parameters of the binary of massive 
black holes. The disadvantage of this or any direct A^-body 
code is the required computational time. However, our calcu- 
lations are accelerated thanks to the special-purpose hardware 
GRAPE-6A single PCI cards of the AEI cluster TUFFSTEIN 
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TABLE 1 

Initial conditions for our featured direct-summation 
w-body models. a/"* is the number of stellar particles used. 
The mass of the binary, normalized to solar masses, is given 

in the THIRD COLUMN, po IS THE INITIAL MAS S DENSITY AT A 
DISTANCE OF 0.1 PC, Wo IS THE KING PARAMETER ( IKlNOll96"6h . ALL 

CASES ARE SINGLE-MASS BUT FOR MODEL C AND F, IN WHICH WE 
HAVE A MASS FUNCTION, SPECIFICALLY A 5-MYRS EVOLVE D KROUPA 
IMF O F MASSES 0.2, 0.5, 50 AND EXPONENTS 1.3 AND 2.3 i IKROUPAI 

120011) . The 7th column shows the initial semi-maior axis of 
the binary in pc. 

used for the si mulations. Each card has a peak performance 
of 130 Gflops (iFukushige et al.ll2005h . so that a single node is 
comparable to a cluster of ~100 individual CPUs working in 
parallel. 

Table [T] gives the initial conditions for the different sim- 
ulations that we feature. We ran six cases with varying 
number densit ies and concentrations, of which two had a 
iKroupal dMTh mass function instead of single-mass stars. 
The IMBHs have equal mass except in simulation F, which 
has a mass ratio of 5. In our simulations the individual time 
steps led to fractional energy errors that were always less than 
10""* per N-body unit of time and, globally, the total energy er- 
ror of the cluster (i.e. the accumulated error in the integration 
of all particles) is 0.015% in the case of our fiducial model 
(A). 

2.2. Evolution of the binary: gravitational radiation versus 
dynamics 

Our approach is to evolve the cluster up to ~ 50-70 Myrs 
using the direct N BODY A code with a central bound binary. 
As we discuss below, we can see in Fig. ([Til that the stel- 
lar cluster experiences a very moderate expansion during the 
hardening of the IMBHs. Once the IMBHs are hard enough 
relative to each other they can be treated as an isolated binary 
that interacts occasionally with a passing star. We note that 
although the subsequent evolution of the binary will not be 
identical to that from the N-body runs, due to the stochastic 
nature of the encounters, the general development is similar 
Such interactions tend to increase the binding energy of the 
binary, hence shrinking its semi-major axis. The eccentricity 
is also changed, both by Newtonian three-body interactions 
fwhi ch can increase or decrease the eccentricitv: see Sesana 
et al. l2007l for a recent treatment) and by gravitational radi- 
ation, which circularizes the binary. The combination of the 
two determines the eccentricity of the binary when it enters 
the frequency range of LISA. 

In Figure |2] we show the inspiral of the binary for mod- 
els A and C. The irregular lines correspond to the A^-body 
simulations. We take the last point of these evolutions and 
the number density of field stars as input to relativistic scat- 
tering experiments_^hich_^e_gerform^^ 
et al. (2006). The equations of motion we use for the three- 
body encounters include relativistic precession to first post- 
Newtonian order, as well as radiation reaction caused by grav- 
itational waves. Between encounters, we evolve the semima- 
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Fig. 1 . — Lagrangian radii showing the evolution of different mass frac- 
tions in the cluster for Model D. The fractions are, from the bottom to the top, 
0.01%, 0.03%, 0.1%>, 0.2%, 0.3% . . . 0.9%) and 0.95% 
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Fig. 2. — Left panel: Inspiral of the IMBH binary of Model E followed 
in the eccentricity-semi-major axis plane. The irregular line shows the re- 
sults of the A?— body simulation. The smooth black solid curves are the esti- 
mated tr ajectories due to gravitational wave emission following the approxi- 
ination of lPetersI 119641) . and the dashed curves show the coiTesponding inspi- 
ral timescale, /qw- The dark dashed area depicts the region of unstable orbits. 
The lightly shaded area corresponds to the phase in the evolution in which 
the ;i = 2 harmonic of the gravitational wave signal is in the LISA band. The 
dashed irregulai" line stalling after the last point of the results of the W-body 
simulation (in the color version depicted in magenta), are the results from 
the scattering expeiiments. See text for further details. Near the beginning 
the eccentricity temporarily exceeds unity because the black holes are not yet 
bound to each other Right panel: Same for Model F, which is one case in 
which we have initially a Kroupa IMF (see Table[TJ. When star with a big 
mass interacts with the binary, the eccentricity change is substantial. The ec- 
centricity therefore wanders up and down, and when it becomes large enough 
the binary has a greater chance to spiral together by gravitational radiation 

jor axis and eccentricity o f the IMBH b inary using the Pe- 
ters quadrupolar formulae (lPeters|[T964l) . The stars that in- 
teract with the binary are sent with a velocity at infinity of 
v=10 km s"', typical of cluster velocity dispersions. The in- 
teraction time is drawn from an exponential distribution with 
a mean time of r = (nEv)"', where « is the stellar number den- 
sity (taken from the A^-body simulations) and S is the scatter- 
ing cross section including gravitational focusing. The typical 
region in which the IMBH binary wanders is larger than its 
radius of influence, hence there is no loss cone as there is for 
supermassive black holes. 
We ran four sets of 40 simulations, two sets starting at large 



Gravitational waves from binary IMBHs 
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Eccentricity 

Fig. 3. — Left panel: Average evolution of eccentricity as a function of 
semimajor axis. Tliis figure sliows the results of 40 three-body scattering 
experiments, starting with IMBH masses of lO' Mq each at an initial semi- 
major axis of 400 AU and an initial eccentricity of 0.55, con'esponding to the 
end point of Model B. We see the -Icr to +lo- ranges of eccentricity. The hor- 
izontal line shows a semimajor axis of 0.093 AU, which is where the orbital 
frequency is 5 X 10"^ Hz and thus the dominant gravitational wave frequency 
is le-4 Hz. 

Right panel: Similar but assuming a Ki'oupa IMF and using the end point of 
Model C. There were no ejections of the binary from the cluster, which was 
assumed to have an escape speed of 50 km/s. The horizontal line is lower 
than in the left panel because the binary mass is less. 

separations with zero eccentricity that established agreement 
with the direct-summation A^-body simulations and two sets 
at the endpoints of Models A and B. In none of these runs was 
the binary itself ejected from the cluster, as expected given 
its large mass. For the Model A endpoint run (binary mass 
600 Mq, initial semimajor axis oo = 10 AU, and initial ec- 
centricity Co = 0.6), the eccentricity was clisa = 0.30 ±0.10 
when the gravitational wave frequency (equal to twice the or- 
bital frequency) was at the lO"** Hz low end of the LISA band; 
for the Model B endpoint run (binary mass 2000 Mq, initial 
semimajor axis aq = 400 AU, and initial eccentricity eo = 0.55) 
we found clisa = 0.24 ± 0.05. We show the envelope of the 
Model B endpoint runs in the left panel of Figure |3] We also 
did scattering experiments corresponding to the endpoint of 
Model C, which had a Kroupa mass function. The results 
are shown in the right panel of Figure |3] Compared to the 
single-mass runs we see considerably greater variance in the 
eccentricity as a function of semimajor axis, and although the 
range of eccentricities in the LISA band clisa = 0.11 ±0.11 
overlaps those in the single-mass runs there are also a num- 
ber of cases in which the binary is nearly circular by the time 
the gravitational wave frequency reaches 10"^ Hz. This could 
be a general feature of scattering interactions when there is a 
broad mass function, but we have not performed enough runs 
to determine this with confidence. 

3. DISCUSSION AND CONCLUSIONS 

We have addressed the inspiral of two massive black holes 
in a single young stellar cluster. Our three main results are: 
(1) the cluster itself experiences only mild structural changes 
as a result of the inspiral, (2) the coalescence takes a short 
enough time (typically < 100 Myrs) that mergers occurs close 
to the time of formation of the cluster, and (3) there is a sig- 
nificant residual eccentricity by the time the binary enters the 
LISA band. We now discuss these conclusions in order. 

The stability of clusters against IMBH mergers is consis- 
tent with analytic expectations even though the binding en- 
ergy of the IMBH binary can exceed the total binding en- 
ergy of the cluster by a significant factor. To see this, con- 
sider a circular IMBH binary of component masses m\ and 
1112 < with total mass M12 = m\+ m -i and r educed mass 
^, = m\m2/M\2. As shown by lOuinlanI (Il996h . a star with 



low speed at infinity that interacts with the binary will typi- 
cally be ejected with a speed Vej ~ . 85 •\/m2 /Mi 2 Vorb , where 
Vorb is the relative speed of the two objects. For equal masses 
m\ = 1112, this is Vq ~ 0.6ybin. Suppose now that the cluster 
has an escape speed Vesc- If Vej < Vesc then the star will be 
retained and share its kinetic energy with the cluster Other- 
wise, the star will be ejected from the cluster without deposit- 
ing significant energy, because the dynamical time of escape 
is much less than the relaxation time (which is the time re- 
quired for the star to give up energy). The binding energy 
of the IMBH binary when v'ej = Vesc will be Eun = ^fJ-Vorh ~ 
5M(^'esc/0.85)^(Mi2/?M2) ~ 0.7miVgs(,. In comparison, if the 
cluster has a three-dimensional velocity dispersion 0-30 and a 
mass Mci, the binding energy of the cluster is E^i ~ jMcictJq. 

The ratio is then fibin/^ci « (mi/Ma) (vesc/o-so)^- Typically 
Ve.sc ~ 2-3 X (T3D, SO Only if the larger black hole mass is 
nil > 0. 1 -0.2Mci could the release of energy unbind the clus- 
ter. We also note that subsequent to this point, the loss of 
mass from stars being thrown out would also soften the clus- 
ter. However, since typically interaction with of order the bi- 
nary mass changes the semimajor axis by a factor of ^ 2, just 
~ 10 Ml 2 in stars will shrink the binary by enough of a fac- 
tor to produce coalescence. Therefore, as verified by our nu- 
merical simulations, hardening of an IMBH binary has only a 
minor effect on the cluster. 

For the time to merger, we note that hardening from large 
separations to a few hundred AU takes ~ 50 Myr, based on 
our simulations. Our three-body runs then indicate that the 
total time from that point to merger is virtually always less 
than 10 Myr, meaning that conservatively the total time from 
formation to merger is less than 10** yr. This is significantly 
shorter than the age of the universe. One consequence of this 
is that if star formation in massive clusters was more common 
at redshift z ^ I than it is now, and if binary IMBH forma- 
tion was also thus more common, then LISA observations of 
IMBH mergers will serve as a unique snapsho t of star forma- 
tion a s well as of cluster dynamics (see also iFregeau et al.l 
l2006h . 

Figure ^ shows that the eccentricity of the binary will be 
in the range ~ 0. 1 - 0.3 when the dominan t gravitat i onal w ave 
frequency is of lO"** Hz. Consistent with lOuinlanI (1 1 996h we 
find that the eccentricity does not undergo a random walk, 
but instead tends to higher eccentricities when the binary is 
hard but before gravitational radiation circularization is im- 
por tant. As discussed in section 4 of Amaro-Seoane & Fre- 
itag ilOO d). a residual eccentricity will induce a difference 
in the phase evolution of the second harmonic compared to a 
circular o rbit, even if it is as small as 0.07. as Amaro-Seoane 
& Freitag (I2006h found. In our case, if we use an eccentric- 
ity ein -4H7 = 0.3 in equation (4) of lAmaro-Seoane & FreitagI 
^2000), we find that the accumulated phase shift A^e > 2tt if 
observations cover a time of at least 

Vg- 0.012 ■(l+z)2yr (1) 

before merger, where z is the redshift. If we set z = 1, then 
we have to cover a time t,^-g = 17 days before merger. This 
means that if we are able to observe the system during that 
period of time before the final coalescence, we will recover 
enough information to determine that the orbit is not circular. 
O n the other hand, if we use a resid ual eccentricity of 0.07, as 
in lAmaro-Seoane & Freitad (l2006l) . we would need 3-4 years 
of observation for a 300 + 300 Mq binary before merger. 
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In conclusion, if young massive clusters form binary 
IMBHs then they will be strong and moderately eccentric 
LISA sources that could serve as unique signposts of clus- 
tered star formation. The non-zero residual eccentricity has 
an impact on the detection of such sources, since it is gener- 
ally assumed that an equal-mass massive binary will have a 
zero eccentricity when entering the LISA band. Our results 
show that in our scenario e is non-negligible for certain cases 
-though for some other models it is very low but detectable-; 
notably, case C and F, which are the only models in which we 
have a mass fraction and thus, they are the more realistic ones. 
The process of formation must of course be studied carefully 
from the standpoints of stellar dynamics and merger product 
evolution, but if binary IMBHs can form then their mergers 
are promising sources for future LISA detections. 
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